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Abstract 

Cluster algorithms are developed for simulating quantum spin systems like 
the one- and two-dimensional Heisenberg ferro- and anti-ferromagnets. The 
corresponding two- and three-dimensional classical spin models with four-spin 
couplings are maped to blockspin models with two-blockspin interactions. 
Clusters of blockspins are updated collectively. The efficiency of the method 
is investigated in detail for one-dimensional spin chains. Then in most cases 
the new algorithms solve the problems of slowing down from which standard 
algorithms are suffering. 



Two-dimensional quantum spin systems are relevant for the description of the un- 
doped anti-ferromagnetic precursor insulators of high-T c superconductors. Presum- 
ably understanding the physics of the superconductors requires an understanding of 
their precursor insulators. This already is a nontrivial problem, which most likely 
does not allow for a complete analytic solution. Therefore it is natural to use a nu- 
merical approach to compute the properties of these materials. At present different 
methods are used in numerical studies of quantum spin systems (for a recent review 
see for example ref.JT|). Small systems can be solved completely by a direct diag- 
onalization of the hamiltonian. For larger systems one can use Monte-Carlo meth- 
ods. For this purpose the finite temperature partition function of the (^-dimensional 
quantum spin system is expressed as a pathintegral of a (d+ l)-dimensional classical 
system of Ising-like spin variables with four-spin couplings. The classical system is 
then simulated on a euclidean time lattice with lattice spacing e using importance 
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sampling techniques like the Metropolis algorithm. This approach was pioneered 
by Suzuki and collaborators 0. An application to one-dimensional spin chains 
is described in ref. |§. For a simulation of the two-dimensional Heisenberg anti- 
ferromagnet see for example ref. ||]. The standard numerical methods are, however, 
slowed down in different ways. First of all most local changes of a configuration 
are forbidden, because many configurations have zero Boltzmann weight. Secondly, 
a change of the magnetization requires a nonlocal update, which has a very small 
acceptance rate at low temperatures. Finally, local algorithms are critically slowed 
down in the continuum limit e — > 0, because then any correlation length in euclidean 
time diverges in units of e. The slowing down of an algorithm is characterized by 
its dynamical critical exponent z. The exponent describes how the Monte-Carlo 
autocorrelation time r oc l/e z (i.e. the time needed to create a new statistically 
independent spin configuration) behaves in the continuum limit. 

For classical spin systems with two-spin couplings — like the Ising model - 
critical slowing down is almost entirely eliminated by the Swendsen-Wang || and 
Wolff || cluster algorithms for which z m 0. Cluster algorithms are, however, 
not directly applicable to models with four-spin couplings. In this paper we map 
the classical spin models with four-spin couplings to blockspin models with two- 
blockspin interactions, which are then simulated using the Swendsen-Wang or Wolff 
cluster algorithms. The blockspin cluster algorithm automatically creates allowed 
spin configurations only. Furthermore, it updates the magnetization efficiently even 
at low temperatures. The efficiency of the blockspin cluster algorithm is investigated 
in detail for the one-dimensional spin chains. We study the critical slowing down by 
investigating the e-dependence of the autocorrelation times of different observables. 
For example, the dynamical exponent of critical slowing down of the algorithm for 
the one-dimensional Heisenberg anti-ferromagnet is z = 0.0(1). The use of improved 
estimators substantially reduces the variance of measured observables. We should 
mention that it is presently not clear how efficient the algorithm works for two- 
dimensional spin systems. 

We consider quantum systems of spins 1/2 with a Hamilton operator 

H = J);S X ■ S x+ jx, (1) 

where S x = a x /2 is a spin operator located at the point a; of a one-dimensional or two- 
dimensional quadratic lattice of even length L with periodic boundary conditions 
and a x are the Pauli matrices. The coupling is between nearest neighbors (fi is the 
unit vector in //-direction) and J is the exchange coupling. For J < parallel spins 
are energetically favored and we have a ferromagnet, while J > corresponds to an 
anti-ferromagnet . 

To express the partition function as a pathintegral the hamiltonian of the one- 



2 



dimensional spin chain is decomposed into H = Hi + H 2 with 



Hi = J Y s x- S x+i , H 2 = J Y S *' S x+i- ( 2 ) 

x=2m X =2m+1 

Now one uses the Trotter formula for the partition function 

Z = Tr exp(-PH) = lim 7V ^ 00 Tr[exp(-e/3i/i) exp(-e(3H 2 )f , (3) 

where e = 1/N is the lattice spacing in the euclidean time direction and /3 is the 
inverse temperature. The original model is recovered in the continuum limit e — > 0. 
For the two-dimensional system one decomposes H — Hi + H 2 + H 3 + H 4 with 

Hi — J X ^ ' S x+ i, H 2 = J Y S x - S x+ i, 

x=(2m,n) x=(2m+l,n) 

H 3 = J S x- S x+2 , H A = J S x- S x+2 , (4) 

x=(m,2n) x=(m,2n+l) 

and one uses the Suzuki formula 

Z = lim iV ^ 00 Tr[exp(-e/3 J ffi) exp(-e[3H 2 ) exp(-e(3H 3 ) exp(-e(3H 4 )] N . (5) 

In the next step the partition function is expressed as a pathintegral of Ising-like 
variables by inserting complete sets of eigenstates |1) and |-1) of a x between the 
factors exp(— ef3Hi). This gives rise to a (d + l)-dimensional classical spin system 
with periodic boundary conditions in the euclidean time direction. The Hi are sums 
of commuting terms. The nonzero elements of the transfer matrix are given by 

(1 1| exp(-e(3JS x ■ SU A )|1 1) = (-1 -1| exp(-e(3JS x ■ S x+fl )\-l -1) = 
exp(-e/3J/4), 

(1 -1| exp(-ef3JS x ■ S x+fl )\l -1} = (-1 1| exp(-ePJS x ■ S x+ii )\-l 1) = 
exp(-e/?J/4)i(l + exp( e /?J)), 

(1 -1| exp(-e(3JS x ■ S x+P )\-1 1) = (-1 1| cxp(-e(3JS x ■ S x+P )\l -1) = 
exp(-e/3J/4)i(l-exp(e/3J)). (6) 

All other elements are equal to zero. For a ferromagnet (J < 0) the expressions are 
positive and can hence be interpreted as Boltzmann factors 

(sis 2 | exp(-ePJS x ■ S x+fl )\s 3 s 4 ) = exp(-S[s u s 2 , s 3 , s 4 ]) (7) 

of spin configurations s = ±1 with a classical euclidean action S. For an anti- 
ferro magnet, on the other hand, this is not possible because some transfer matrix 
elements are negative. Therefore in the anti-ferromagnetic case one performs a 
unitary transformation of the original hamiltonian by rotating every second spin by 
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an angle it. The nonzero elements of the transfer matrix are then positive for J > 0. 
Expressed as a pathintegral the partition function takes the form 



Z = U £ exp(-S) 

x,t s(x,t)=±l 



with 



x 



x 



exp(— S) = j~[ exp(—S[s(x,t),s(x + l,t),s(x,t + l),s(x + i,t + l)]) 

x=2m,t=2p 

x J| exp(—S[s(x,t),s(x + i,t),s(x,t + l),s(x + i,t + l)]) 

x=2m+l,t=2p+l 

(9) 

in the one-dimensional case and with 

exp(— S) = Yl exp(—S[s(x,t),s(x + l,t),s(x,t + l),s(x + l,t+l)\) 

x=(2m,n),t=Ap 

Yl exp(-S[s(x, t), s(x + i, t), s(x, £ + 1), s(x + i, t + 1)]) 

x=(2m+l,n),t=4p+l 

IJ exp(-5[s(x, £), s(x + 2,t), s(x, t + l),s(x + 2,t+ 1)]) 

x=(m,2n),t=4p+2 

x Yl exp(-S[s(x,t),s(x + 2,t),s(x,t + l),s(x + 2,t+l)\) 

x=(m,2n+l),t=4p+3 

(10) 

for the two-dimensional system. Note that the classical spins s(x,t) interact with 
each other via four-spin couplings. Most spin configurations are forbidden (their 
Boltzmann factor vanishes) because the corresponding elements of the transfer ma- 
trix are zero. This is a problem for standard local algorithms because most local 
changes of a configuration are not allowed. Therefore it is natural to attempt a 
collective nonlocal update of the spins. For spin models with two-spin couplings 
this can be done using the Swendsen-Wang or Wolff cluster algorithms, which 
flip whole clusters of spins simultaneously. The cluster algorithms can, however, not 
be applied directly to models with four-spin interactions. 

To make an application of cluster algorithms possible we map the classical spin 
models with four-spin couplings to blockspin models with two-blockspin couplings. 
A blockspin is a collection of a few spins. For the one-dimensional spin chain we 
define blockspins consisting of four spins located at the corners of a plaquette 

6(2m, 2p) = {s(x, t), s(x + l,t), s(x, t + 1), s(x + 1, t + 1)} (11) 

with x = 2m— 1 and t = 2p. Note that each spin belongs to exactly one blockspin and 
the blockspins live one a lattice with a doubled lattice spacing. The original action 
determines the action for the blockspins. For example, the four-spin interactions for 
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Figure 1: Three different blocking schemes. The blockspins b and b consist of four 
spins and live on lattices with a doubled lattice spacing. The crossed plaquettes 
carry the four-spin interaction which turns into a two-blockspin interaction. In the 
third scheme one blockspin b consists of all spins with the same space coordinate 

Xq. 




x = 2m, t = 2p induce couplings between space-like nearest neighbor blockspins at 
(2m, 2p) and (2m + 2, 2p) 

S[b(2m, 2p), b{2m + 2, 2p)\ = S[s(x, t), s(x + 1, t), s(x, t + 1), s(x + i, t + 1)], (12) 

and the four-spin interactions for x = 2m + 1, t — 2p + 1 induce couplings between 
time-like nearest neighbor blockspins at (2m, 2p) and (2m, 2p + 2) 

S[b{2m,2p),b{2m,2p + 2)] = S[s(x, t), s(x + i, t), s(x,t + 1), s(x + 1, t + 1)]. (13) 

The spins can also be arranged to blockspins in another way 

b(2m, 2p) = {s(x, t), s(x + 1, t), s(x, t + 1), s(x + i, * + 1)} (14) 

where now x = 2m and t = 2p — 1. Then the four-spin interactions for x = 
2m, t = 2p induce two-blockspin couplings between time-like nearest neighbors and 
the interactions at x — 2m + l,t = 2p + 1 induce couplings between space-like 
nearest neighbors. The two blocking schemes of spins into blockspins are illustrated 
in figs.la,b. To ensure ergodicity an updating algorithm for the blockspins must 
alternate between the two blocking schemes. 

For the two-dimensional system a blockspin consists of eight spins located at the 
corners of a cube 

b(2m, 2n, 2p) = {s(x, t),s(x + l,t), s(x + 2, t), s(x + 1 + 2, t), 

s(x,t + l),s(x + i,t + l),s(x + 2,t + l),s(x + i + 2,t + l)} (15) 

with x = (2m — l,2n — 1) and t = 2p. Again the blockspin lattice has a doubled 
lattice spacing. Another blocking scheme is given by 

b(2m, 2n, 2p) = {s(x, t),s(x + l,t), s(x + 2, *), s(x + 1 + 2, t), 

s(x,t + l),s(x + i,t + l),s(x + 2,t + l),s(x + l + 2,t + l)} (16) 
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where now x = (2m, 2n) and t — 2p — 1. Using the action of eq.([T0|) it is straight- 
forward to show that both schemes induce two-blockspin interactions only. Again, 
to ensure ergodicity the algorithm should alternate between the two schemes. 

The cluster algorithms make use of a flip symmetry of the blockspin model. A 
blockspin b = {s 1; s 2 , s n } is flipped to —b = {— si, — s 2 , — s n } simply by flipping 
all spins that belong to b. It is clear that the blockspin action is invariant against 
flipping all blockspins simply because the original action is invariant against flipping 
all spins. The cluster algorithm puts bonds between nearest neighbor blockspins b 
and b' with the probability 

p= 1 -min{l,exp(-5[-6,6 / ])/exp(-5[6,6 / ])}. (17) 

Two blockspins which are connected by a bond belong to the same cluster. The 
algorithm flips all blockspins in one cluster simultaneously. In the Swendsen-Wang 
multi-cluster method each cluster is flipped with the probability 1/2. In the Wolff 
single-cluster method one blockspin is randomly selected and the cluster to which it 
belongs is flipped with probability 1. Recently, it has been realized that the single- 
cluster algorithm can be vectorized 0. Using S[b, —b'] = S[—b, b'] and S[—b, —b'] = 
S[b, b'\ one can show that both algorithms obey detailed balance. The blockspin 
cluster algorithm automatically generates allowed spin configurations only. If one 
attempts to flip b when exp(— S[— b, b'}) = 0, i.e. when the new configuration is 
forbidden, the algorithm puts a bond with probability p = 1 — min{l, 0} = 1. Hence 
b and b' fall in the same cluster and must be flipped simultaneously. The cluster 
growth comes to an end only when the new configuration is allowed. When one 
uses the single-cluster method it is guaranteed that one generates a new different 
configuration in each step. 

As defined up to now the blockspin cluster algorithms are not ergodic because 
they change the magnetization M — 1/2 J2 X s i x i to) and the staggered magnetization 
M s = 1/2 Ei( _ l)M x ^o) by even numbers only. M and M s are integers (rather 
than half-integers) because we work on a lattice with an even number of space points. 
Note that M is independent of to because the transfer matrix commutes with the 
total spin. To allow M and M s to change by odd numbers also, we introduce extra 
blocking schemes with some block spins bo = {s(xo,t) for all t} consisting of all 
spins with the same space coordinates xo- The other blockspins are taken from the 
above schemes. A typical example is shown in fig.lc. Similarily, we include blocking 
schemes which allow the so-called winding number N w = 1/2 Ya s ( x 0i to change 
by an odd integer. 

The Swendsen-Wang and Wolff cluster algorithms eliminate critical slowing down 
from the Ising model because it is not frustrated. For strongly frustrated models, 
for example for spin glasses, cluster algorithms do not work efficiently Therefore 
the question arises if our blockspin models are frustrated or not. To answer this 
question for the spin chains one must investigate all allowed configurations of four 
nearest neighbor blockspins located at the corners of a square of the blockspin lattice. 
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Nearest neighbor pairs of blockspins interact via bonds (the four-spin interactions 
of the original classical spin model). Now the Boltzmann factor of a bond is com- 
pared to the value it has when one of the two blockspins at its ends is flipped. 
The blockspin configuration with the larger Boltzmann factor is the favored config- 
uration of this bond. If one can flip the four blockspins such that all four bonds 
between them are in their favored configuration the configuration is not frustrated. 
We have verified that the blockspin model for the one-dimensional anti-ferromagnet 
has indeed no frustrated configurations. For the one-dimensional ferromagnet, on 
the other hand, some configurations are frustrated. In a frustrated configuration 
one of the four bonds cannot be in its favored configuration. However, one can 
show that the Boltzmann factor of the disfavored bond is in all cases only a factor 
(l + exp(e/3J))/2 smaller than the favored Boltzmann factor. In the continuum limit 
e — > the factor goes to 1 and the frustration disappears. A priori it is not clear if a 
weak frustration can slow down the cluster algorithm in the ferromagnetic case. As 
we will see later, the blockspin cluster algorithm does in fact work perfectly only in 
the anti-ferromagnetic case. However, also for the ferromagnet the autocorrelation 
times of the cluster algorithm are much smaller than the ones of the Metropolis al- 
gorithm. For the two-dimensional spin systems both for the ferromagnet and for the 
anti-ferromagnet some blockspin configurations are frustrated. In addition, not all 
frustrations disappear in the continuum limit. Therefore, one should not expect that 
the blockspin cluster algorithm for two-dimensional quantum spin systems works as 
well as in the one-dimensional case. However, we expect that it still decorrelates 
faster than standard algorithms. 

Cluster algorithms offer the possibility to use improved estimators which reduce 
the variance of different observables. For example, for the single-cluster method the 
susceptibility can be expressed as 

R M 2 

x = f- d ( M2 ) = 2dN ^)> ( 18 ) 

where L d is the spatial volume, 2dN is the number of points in the euclidean time 
direction, \C\ = Z)(x,t)ec 1 is the cluster size and M c = l/2J2(x,t )ec s ( x ^o) is the 
cluster magnetization. It is interesting to note that also Mc is independent of t - As 
a consequence, clusters with nonzero magnetization must wrap around the lattice 
in the euclidean time direction. Small clusters which do not wrap around the lattice 
have Mc = 0. Let L d be the minimal spatial extent of a cluster closed in euclidean 
time, i.e. L d = mm to {J2 ( x ,t )ec !}• The cluster magnetization is then limited by 
\Mc\ < Lc/2 and the cluster size is restricted by \C\ > 2dNL d such that 

(L d c ) > A -f- (19) 

The algorithm updates the susceptibility efficiently only if many clusters wrapping 
around the euclidean time direction fit into the lattice. This requires (L d ) <C L d . 
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Hence we expect that the cluster algorithm works efficiently only if 



X « (20) 

The ferromagnetic systems have degenerate ground states with total spin L d /2. In 
the zero temperature limit the susceptibility is then given by x — fl{L d + 2)/12 and 
the algorithm does not work efficiently. In all other cases, for example at non-zero 
temperature or for anti-ferromagnets, x stays finite as the spatial volume goes to 
infinity and the algorithm should work. One can also define an improved estimator 
for the staggered susceptibility Xs — P(Mg) /L d . Since the staggered magnetization 
is to-dependent also small clusters contribute to Xs- Finally, we introduce the internal 
energy density e = —1/L d (d \nZ/dj3). 

Now let us turn to the numerical results. We have tested the blockspin cluster 
algorithm in detail for one-dimensional spin chains. To test the efficiency of the 
blockspin cluster algorithm we compare it to a Metropolis update of the blockspins. 
In fact, the blockspin Metropolis algorithm is similar to the standard algorithms 
commonly used. For both algorithms we measure the autocorrelations of the sus- 
ceptibility 

C x (5r) = x(r)x(r + 5r) (21) 

as a function of the Monte-Carlo time difference Sr. From this we obtain the inte- 
grated autocorrelation time r x as 

oo oo 

exp(-l/r x ) = £ C x (St)/ C x (5r). (22) 

5t=1 8t=0 

Note that r x = r for C x (5r) oc exp(— 5t/t ). In the same way we define the 
integrated autocorrelation times t Xs of the staggered susceptibility and r e of the 
internal energy density. For the cluster algorithm we use the single-cluster method 
and we measure the autocorrelations based on the improved estimators. Our results 
are summarized in table 1. The quoted autocorrelation times are in units of Monte- 
Carlo sweeps. A sweep is defined such that it takes about the same CPU-time for 
the cluster and for the Metropolis algorithm. In all cases we have performed a 
random start followed by 5000 sweeps for thermalization and by 50000 sweeps for 
measurements. 

The Metropolis algorithm can change the magnetization only by flipping a block- 
spin b which consists of all spins with the same space coordinate x . The flip is 
allowed only if all spins of &o are parallel. Surprisingly, the probability for this 
situation — and hence the acceptance rate of the global step — is finite in the 
continuum limit. Therefore the autocorrelation times r x and t Xs do not diverge in 
the continuum limit. However, they are at least an order of magnitude larger than 
the ones of the cluster algorithm. The autocorrelation time of the internal energy 
density, on the other hand, diverges as r e oc l/e Ze . For the Metropolis algorithm 
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Table 1: Results for the one-dimensional Heisenberg spin chains. The blockspin 
single-cluster algorithm (C) is compared to the Metropolis algorithm (M). 
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our (3 = 1 data both for the ferro- and for the anti-ferromagnet yield a dynamical 
critical exponent z e = 0.8(1) while the cluster algorithm has z e = 0.0(1). Note that 
the exponent may depend on the observable considered, because it was extracted 
from the integrated autocorrelation time and not from the exponential decay of an 
autocorrelation function. 

Although the global step of the Metropolis algorithm is not critically slowed 
down (its acceptance rate is e-independent) it is very severely slowed down at low 
temperatures (large 0). This prevents the application of standard numerical meth- 
ods to quantum spin systems at low temperatures. Already at (5 = 8 the Metropolis 
algorithm has r x = 3300(200) for the ferromagnet and r x = 200(20) for the anti- 
ferromagnet. The cluster algorithm, on the other hand, has autocorrelation times of 
at most a few sweeps and one can go to lower temperatures like /? = 16 very easily. 
For the ferromagnet the inequality x ^ (3L d /A puts a lower limit on the temperature 
when the spatial volume is fixed. Our low temperature data for the ferromagnet at 
/? = 16 show some slowing down but the autocorrelation times are moderate. For 
the anti-ferromagnet, on the other hand, there is no indication of slowing down 
when the temperature is lowered or when the continuum limit is approached and 

Z X = Z Xs = Z e = 0.0(1). 

To summarize we have developed blockspin cluster algorithms for one- and two- 
dimensional quantum spin systems. Our results for spin chains show that the new 
algorithms eliminate slowing down for anti-ferromagnets. For ferromagnets slowing 
down is also practically eliminated as long as the temperature is not too small. One 
should, however, be careful in generalizing our findings to two-dimensional quan- 
tum spin systems. Since the corresponding blockspin models have some frustrated 
configurations, it is presently not clear how efficient the algorithm works in this 
case. A detailed analysis of the efficiency of the blockspin cluster algorithm for 
the two-dimensional Heisenberg anti-ferromagnet is in progress. We like to mention 
that blockspin cluster algorithms can also be applied to other models. First of all 
the algorithm also works for spin systems with anisotropic couplings or (with slight 
modifications) in an external magnetic field. Secondly, it can be applied to models of 
higher spins (S = 1, 3/2, ...) if it is combined with a method that changes the abso- 
lute value of the spin projection \S^\. The one-dimensional Heisenberg spin chain is 
equivalent to a six- vertex model which is a special eight- vertex model. General eight- 
vertex models with a spin flip symmetry can also be treated with blockspin cluster 
algorithms. Finally, the algorithm may be useful for simulations of one-dimensional 
Fermi-systems like the ones described in ref. ||. 

It is a pleasure to thank P. Hasenfratz who initiated our interest in quantum 
spin systems for his support. 
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